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ABSTRACT 


Several  ocean  circulation  models  are  evaluated  and  compared  on  coastal  applications. 
The  models  are  all  primitive  equation,  tree-surface  models,  that  offer  a  wide  range  of 
mathematicai  and  numerical  formulations. 

Two  potential  problem  areas  have  previously  been  identified  with  the  type  of  the  com¬ 
pared  ocean  models:  the  use  of  cr-coordinates  in  steep  bathymetry  and  unrealistically 
high  dissipation  that  may  prevent  the  formation  of  high  gradient  regions,  simulation  of 
small-scale  features,  and  the  proper  onset  of  dynamic  instabilities.  A  principal  focus  of 
this  research  is  to  evaluate  the  degree  to  which  these  problems  limit  the  performance  of 
the  models  in  coastal  applications.  The  model  comparison  is  conducted  over  a  series 
of  experiments  on  semi-idealized  configurations. 

This  study  indicates  that  the  solutions  may  sensibly  differ  from  model  to  model.  This 
is  primarily  due  to:  1)  different  response  to  the  calibration  of  the  physicai  and  numeri¬ 
cal  parameters,  2)  perturbations  in  the  dynamical  balances.  In  general,  differences  in 
the  primary  features  (such  as  front  location,  strength  of  the  coastal  jet)  can  be  adjusted 
by  adequate  calibrations.  Representation  of  the  secondary  features  (such  as  internal 
wave  propagation,  temperature  and  salinity  anomalies)  may  depend  upon  different  pa- 
rameterizations  and  numerical  treatments  of  the  same  physical  terms  that  may  alter  the 
relative  importance  of  the  dominant  processes. 
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GOALS 


The  long-term  goal  is  to  contribute  to  the  development  of  ocean  predictive  systems  to 
support  Naval  operations  in  coastai  areas.  To  date,  the  predictive  capability  in  shallow 
water  has  not  been  addressed  vigorously,  and  several  aspects  are  still  unresolved.  This 
research  will  support  the  choice  of  the  ’most  suitable  and  accurate  model  formulations 
to  be  applied  in  operational  applications  over  the  littoral  areas. 


OBJECTIVES 

The  objectives  are  to  evaluate  several  ocean  circulation  models  configured  on  limited- 
area  coastal  domains,  compare  the  physical  and  numerical  formulations,  and  prepare 
a  series  of  guidelines  and  recommendations. 


IMPACT 

With  the  shift  of  the  Navy’s  environmental  interests  from  deep  water  to  coastal  regions, 
there  is  the  need  to  redirect  modeling  requirements.  Quantities  of  importance  to  Naval 
coastal  operations  include  water  depth,  current  speed  and  direction,  surf,  waves  and 
swell,  thermal  and  salinity  structure,  and  water  clarity.  This  study  primarily  focuses  on 
the  prediction  of  water  depth,  current  speed  and  direction,  and  the  evolution  of  the  ther¬ 
mal  and  salinity  fields. 
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1.  INTRODUCTION 


In  past  years,  the  Navy’s  primary  ocean  modeling  efforts  have  been  focused  on  devel¬ 
oping  global  or  basin-wide  3-D,  eddy-resolving  ocean  circulation  models.  With  the  shift 
of  the  Navy’s  environmental  interests  from  deep  water  to  coastal  regions,  there  is  the 
need  to  redirect  modeling  requirements.  It  is  recognized  that  coastal  dynamics  is  physi¬ 
cally  quite  different  from  the  deep  ocean.  For  example,  typical  time  scales  over  the  shelf 
are  generally  much  shorter  than  on  deep  waters.  Thus,  prediction  on  coastal  regions 
requires  a  different  approach  than  forecasting  at  basin  scale. 

It  is  important  to  separate  the  physical  regimes  that  can  be  realistically  modeled,  from 
those  for  which  the  model  formulations  are  unsuitable.  Test  experiments  are  needed  to 
compare  model  performances  and  evaluate  the  accuracy  of  the  solutions.  Quantities  of 
importance  to  Naval  coastal  operations  include  water  depth,  current  speed  and  direction, 
surf,  waves  and  swell,  thermal  and  salinity  structure,  and  water  clarity.  This  research 
will  primarly  focus  on  the  prediction  of  water  depth,  current  speed  and  direction,  and  the 
evolution  of  the  thermal  and  salinity  fields. 

The  approach  is  to  compare  existing  ocean  models  in  limited-area,  coastal  applications. 
Four  free-surface,  primitive  equation  models  have  been  selected  for  the  comparison: 

•  the  Princeton  Ocean  Model  (POM,  Mellor  and  Blumberg), 

•  the  Sigma  Z-level  Model  (SZM,  Martin), 

•  the  S-Coordinate  Rutgers  University  Model  (SCRUM,  Song  and  Haidvogel), 

•  the  Estaurine  Coastal  Ocean  Model  Semi-Implicit  (ECOMSI)  (Blumberg). 

The  models  offer  a  wide  range  of  physical  formulations  (such  as  lateral  and  vertical  mix¬ 
ing),  and  numerical  algorithms  (such  as  spatial  and  temporal  schemes,  implicit/explicit 
treatment  of  the  barotropic  flow)  (Tab  1 .). 

Model  comparisons  offer  a  number  of  benefits  and,  to  some  degree,  support  the  choice 
^  of  the  ’’most  suitable  and  accurate”  model  to  be  transferred  to  operational  applications, 
while  shedingTight  on  its  forecasting  skills.  The  comparisons  can  demonstrate  the  su¬ 
periority  of  particular  model  formulations,  model  parameterizations,  and  numerical  tech¬ 
niques.  The  scatter  of  the  results  from  models  that  have  been  determined  to  be  ade¬ 
quately  and  correctly  formulated  may  provide  an  indication  of  the  degree  of  uncertainity 
to  be  expected  in  the  model  predictions.  Features  that  have  been  found  to  be  partic¬ 
ularly  efficient  and/or  accurate  can  be  extracted  from  one  model  and  incorporated  in 
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another. 


If  two  model  solutions  are  significantly  different,  the  differences  can  be  examined  to 
determine  why  the  models  differ,  and  which  model  Is  most  likely  to  be  correct.  If  mod¬ 
els  agree,  it  does  not,  guarantee  that  they  are  both  correct,  but  it  does  increase  our 
confidence  in  the  results. 

The  physical  accuracy  and  the  numerical  character  of  the  models  are  evaluated  on  a 
variety  of  test  problems.  The  coastal  study  processes  include  1)  wave  propagation,  2) 
upwelling  and  downwelling  regimes,  3)  set  up  of  wind  driven  circulation,  and  4)  genesis 
and  evolution  of  fronts  instabilities. 


2.  MODEL  ASSESSMENTS 

The  models  have  generally  been  retrieved  by  anonymous  ftp.  With  the  exception  of 
SCRUM,  the  programs  were  usually  condensed  in  one  or  two  file.  That  made  it  difficult 
to  isolated  parts  of  the  codes  or  verify  parameter  values.  Therefore,  all  the  models 
have  been  configured  on  a  modular  form  and  compiled  through  the  ’make  command. 
Following  SCRUM’S  model  architecture,  C-preprocessor  derectives  were  included  in 
the  F77  programs.  This  allows  one  to  conditionally  include  part  of  the  code  before  the 
F77  compile  sees  them,  to  save  the  expense  of  compiling  and  linking  options  that  are 
not  being  used  in  a  given  model  configuration,  and  to  add/bypass  parts  of  the  code 
with  minimum  changes  in  the  programs.  The  modifications  have  increased  the  models’ 
portability  on  several  platforms,  and  reduced  the  compiling  process  when  only  a  few 
changes  were  made  in  a  model’s  configuration. 

All  the  models  were  also  linked  to  the  plotting  package  developed  by  K.  Hedstrom  at 
Rutgers  University,  which  makes  use  of  the  NCAR  Graphic  software.  This  made  it  pos¬ 
sible  to  visually  compare  model  solutions  on  a  common  ground  for  a  more  stringent 
evaluation.  Since  the  variability  of  the  model’s  solutions  is  one  of  the  major  topics  of  the 
comparison,  the  output  has  been  interfaced  with  the  Empirical  Orthogonal  Functions 
(EOF)  programs  developed  by  D.  Fox  at  NRL  and  recently  implemented  at  COAM. 

The  models  are  configured  with  idealized  geometries  and  with  sufficient  horizontal  and 
vertical  resolution  to  provide  accurate  estimates  of  the  dynamics.  All  the  tests  are  con- 
fgured  on  a  periodic  channel  to  avoid  the  application  of  open  boundary  algorithms  that 
could  have  altered  the  model  evaluation  and  comparison. 
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3.  ROBUSTNESS 


One  issue  addressed  is  to  determine  the  robustness  and  portability  of  the  models.  Sim¬ 
ulations  investigate  how  and  if  mixing  algorithms,  filtering  and  smoothing  procedures 
contribute  to  the  models’  numerical  stability  at  the  expenses  of  the  physical  accuracy. 


3.1  The  Choice  of  the  Vertical  Coordinate 

The  choice  of  the  vertical  coordinate  transformation  is  a  controversial  and  unresolved 
issue.  The  ^-levels  have  the  flexibility  of  a  high  resolution  near  the  surface,  for  a  good 
representation  of  the  Mixed  Layer  (ML)  dynamics,  and  a  coarser  resolution  in  the  deeper 
strata,  so  that  they  are  numerically  efficient.  However,  since  the  bathymetry  is  approx¬ 
imated  by  a  step-like  function,  topographic  effects  may  be  poorly  parameterized  (Mar¬ 
tin,!  997  manuscript  in  preparation). 

A  series  of  problems  have  been  recognized  and  documented  for  the  (j-coordinates.  The 
pressure  forces  and  vertical  advection  terms  may  be  affected  by  the  computer  truncated 
arithmetic,  with  sensible  loss  of  accuracy  and  stability.  For  n  evenly-spaced  a-levels, 
the  hydrostatic  consistency  is  defined  by  the  ratio: 

1  AH 
nH  Ax 

(Mellor  et  al.,1994;  J.  Atm.  and  Ocean.  Tech.).  Although  the  problem  may  become 
severe  in  coarse  grid  resolution  with  steep  topography  or/and  strong  stratification,  the 
topography  gradient  and  the  numerical  grid  resolution  for  most  of  the  coastal  applica¬ 
tions  satisfy  the  hydrostatic  consistency  condition. 


3.2.  The  Implicit  and  Explicit  Formulations 

It  is  well  known  that  Implicit  schemes  are  less  accurate,  but  more  efficient  than  the 
explicit  schemes  (Grotjan  and  O’Brien,  1976;  Mon.  Weath.  Rev.).  Explicit  schemes  are 
always  more  accurate  in  wave  propagation  and  energy  trasport,  even  for  those  waves 
that  are  solved  by  several  grid  points.  Implicit  schemes  retain  stability  by  slowing  down 
the  traveling  features. 

There  is  a  fundamental  difference  between  implicit  and  explicit  treatment  of  the  barotropic 
flow.  Explicit  schemes  have  a  local  character:  at  any  grid  point  the  solution  is  updated 
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by  knowledge  of  the  fields  at  the  previous  time  steps  and  closest  surrounding  points. 

On  the  other  hand,  implicit  schemes  depend  upon  the  application  of  algorithms  for  the 
solution  of  2-D  elliptic  equations,  which  present  a  global  character:  the  field  is  obtained 
by  solving  simoultaneously  an  algebraic  set  of  equations  that  connects  points  geogra- 
phycaliy  distant. 

3.2.1.  The  Truncation  error 

Elliptic  solvers  usually  rely  on  iterative  methods  and  a  compromise  has  always  to  be 
found  between  accuracy  and  computational  cost.  Experiments  evaluate  how  the  barotropic 
flow  is  affected  by  the  choice  of  the  truncation  error  in  an  iterative  procedure.  The  tests 
are  conducted  with  the  SZM  model  which  includes  the  Simuoulaneous  Over-Relaxation 
(SOR)  cyclic  boundary  algorithm.  Although  the  SOR  procedure  is  most  likely  to  not  be 
used  in  realistic  coastal  applications,  it  makes  evident  problems  that  it  may  share  with 
other  more  sophisticated  algorithms.  The  simulations  analyze  the  response  to  a  long¬ 
shore  wind  forcing  of  a  stratified  channel  with  an  uniform  along-channel  bathymetry,  so 
that  the  problem  is  truly  2-D.  The  sea  surface  mean  values  appear  to  be  a  good  indi¬ 
cator  as  to  whether  or  not  the  solver  has  reached  a  reasonable  accuracy.  The  physical 
problem  conserves  the  initial  mass,  but  the  computed  mean  surface  values  are  roughly 
of  the  same  order  of  magnitude  of  the  specified  truncation  error  after  1  day  of  simula¬ 
tions.  Although  differences  in  the  velocity  values  are  within  errors  of  few  millimeters  per 
second,  the  cumulative  truncation  error  may  sensibly  change  the  current  systems  over 
extended  simulations. 


3.2.2.  The  correction  procedure 

The  presence  of  the  nonlinear  terms  in  the  barotropic  flow  equations  implies  terms  such 
as  h“+iu®+^  that  complicate  the  application  of  the  implicit  schemes.  In  general,  the 
advective  terms  are  approximated  by  h“u“+^,  in  order  to  avoid  the  solution  of  a  nonlin- 
ea;  differential  equation.  However,  the  formulation  may  lead  to  inconsistencies  of  the 
numerical  algorithm. 

SZM  introduces  an  iterative  procedure  to  correct  the  advection  terms  of  the  momen¬ 
tum  equations.  It  is  found  that  the  correction  procedure  increases  the  accuracy  of  the 
semi-implicit  formulation.  At  the  second  corrective  iteration,  the  number  of  the  SOR’s 
iterations  may  range  from  20  to  30%  of  the  previous  passage,  indicating  that  at  least 
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two  passages  are  necessary  to  achieve  a  higher  convergence  of  the  elliptic  solver,  and 
a  higher  consistency  of  the  implicit  scheme. 

Although  the  simulations  support  the  usefulness  of  a  corrective  algoritm,  the  procedure 
is  neglected  in  most  of  the  free-surface,  semi-implicit  models.  The  weakness  of  the  SZM 
approach  is  the  computational  cost,  since  the  corrective  iterations  are  performed  over 
the  2-D  and  3-D  fields.  It  would  be  important  to  verify  the  efficiency  and  accuracy  of  a 
corrective  procedure  applied  to  the  external  mode  only. 


3.3.  Time  Filtering 

Both  POM  and  SZM  make  use  of  the  Asselin  filter  to  avoid  the  splitting  modes  of  the 
leapfrog  scheme.  As  Appendix  A  illustrates,  the  Asselin  filter  is  diffusive.  Consider  a 
wave  that  travels  one  grid  spacing  over  fj,  time  iterations,  such  that  C  =  Ax/(^At).  For 
the  vaiues  Ax  =7.5  km.  At  =  900  s,  and  the  typical  Asselin  parameter  a=  0.05,  the 
diffusive  coefficient  of  (A.5)  is:  Nc  ~  1 500/ 

i)  slow  waves :  //  =  10  iVc  ~  9  m^s“^ 

ii)  fast  waves :  /i  =  2  Nc  ~450  m^s"^ 

For  slowly-varying  (with  respect  to  the  time  increment)  features,  time  spiitting  is  not 
a  severe  problem  and  generally,  the  computational  mode  does  not  amplify  to  criticai 
levels  (Mesinger  and  Arakawa,  1976;  GARP  Series,  17),  and  the  Asselin  filter  may 
not  be  an  essential  component  for  the  model  stability.  However,  for  the  fast  waves,  the 
computational  eddy  coefficient  of  the  Asselin  filter  may  became  sensibly  greater  than  the 
average  value  of  the  Smagorinsky  coefficients,  and  the  algorithm  may  alter  the  dominant 
dynamicai  balances  of  the  fast  traveling  features. 


3.4.  Horizontal  Diffusion 

'  One  of  the  common  potential  problems  for  ^-coordinate  models  is  the  application  of  the 
horizontal  diffusion  operator  acting  on  surfaces  of  constant  sigma.  The  numerical  for¬ 
mulation  is  an  approximation  of  the  diffusion  operator  on  the  ^-transformation  and  might 
generate  computational  heating/cooling  sources  in  regions  of  steep  topography.  Exper¬ 
iments  evaluate  the  physical  accuracy  and  computational  efficiency  of  the  geopotential 
diffusion  operator  introduced  by  Barnier  et  al.  (1 996,  Deep  Sea  Res.).  As  Fig.  1 .  illus¬ 
trates  the  rotation  of  the  tensor  eliminates  the  spurious  heating  in  the  proximity  of  the 
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continental  break  and  shelf.  The  new  algorithm  sensibly  reduces  the  across-channel 
gradient  in  the  constituent  fields  at  the  shelf  break. 

3.5.  Hybrid  Numerical  Schemes 

The  early  versions  of  the  SCRUM  model  (henceforth  referred  as  SCRUMJ)  were  de¬ 
signed  with  hybrid  schemes,  such  as  II  order  finite  differencing  for  the  horizontal,  and  I 
order  finite  element  for  the  vertical  derivative  terms.  In  this  configuration,  the  model  did 
not  performe  at  the  expected  level.  The  newest  version  of  the  SCRUM  model  (hence¬ 
forth  referred  as  SCRUM_^)  has  removed  the  vertical  finite-element  formulation  and 
substituted  it  with  the  classical  II  order  finite  differencing,  and  the  models’  preformance 
has  improved  tremendously.  This  may  suggest  that  it  is  important  that  terms  that  are 
dynamically  in  balance  be  numerically  treated  with  schemes  at  least  with  the  same  trun¬ 
cation  errors.  Unfortunately,  SCRUM_2  has  been  released  toward  the  end  of  this  project 
and  has  been  tested  in  only  few  cases. 

4.  MODEL  COMPARISON 
4.1.  Wave  Propagation 

Waves  are  an  essential  component  of  coastal  dynamics  and  provide  an  useful  test  prob¬ 
lem  for  model  evaluation.  The  classic  freely  propagating  linear  waves  are  analyzed: 
barotropic  shelf  waves  (BW),  Kelvin  (KW)  and  Poincare  waves  (PW).  For  those,  an 
analytical  solution  is  known  in  the  rigid  lid  (BW)  and  small  amplitude  (KW  and  PW)  ap¬ 
proximations.  The  model  accuracy  is  determined  by  the  correlation  factor  and  rms  error 
(in  percent)  between  the  numerical  and  analytical  solutions.  The  goals  of  the  experi¬ 
ments  are  twofold:  they  seek  1 )  to  determine  the  accuracy  of  the  solutions  as  a  function 
of  number  of  grid  points  per  wavelength  and  2)  to  analyze  how  advection  and  mixing 
terms  may  distort  the  propagation  of  the  free  waves.  Thus,  the  experiments  are  con¬ 
figured  also  as  a  function  of  the  Rossby  number  e  =  UL/i.  Henceforth,  the  term  linear 
refers  to  simulations  in  which  the  wave  amplitude  has  been  specified  to  keep  the  Rossby 
number  small. 

All  the  models  show  high  correlation  factors  between  the  analytical  and  numerical  solu¬ 
tions  for  the  linear  waves  solved  by  several  grid  points.  The  correlation  factor  decreases 
as  the  number  of  grid  points  per  wavelength  is  reduced.  As  the  number  of  grid  point  per 
wavelenght  is  reduced  or  the  Rossby  number  is  increased,  the  correlation  factors  de- 
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crease  rapidly  for  all  the  models,  but  without  sensible  differences  among  the  models’ 
performance. 


4.2.  Long-term  Wind  Driven  Circulation 

The  tests  analyze  the  response  to  an  along-shore  wind  on  a  periodic  channel.  The 
channel  is  377  km  wide  and  535  km  long.  The  shelf-break  is  parabolically  shaped.  The 
shallow  region  (30  m  depth)  is  30  km  wide  at  the  middle  channel  and  100  km  at  the 
boundaries.  The  maximum  depth  is  300  m  and  the  mean  bathymetry  gradient  is  0.005. 
The  wind  stress  (1 .2  m^s"^)  is  steady  and  uniform  in  the  along-shore  direction  and  set 
to  zero  in  the  deep  part  of  the  channel  with  a  cos-like  function. 


4.2.1  Downwelling 

POM  generates  the  smoothest  and  weakest  coastal  jet.  The  maximum  along-channel 
velocities  are  at  the  shelf  break  and  the  coastal  jet  extends  well  over  the  whole  continen¬ 
tal  slope.  SZM’s  coastal  jet  is  narrower  and  located  more  offshore,  with  the  axis  in  the 
deepest  part  of  the  continental  slope.  The  coastal  current  system  of  SCRUMJ  solution 
is  located  further  offshore,  and  it  is  generally  wider  than  in  the  other  models.  SCRUMJ 
also  produces  the  highest  horizontal  velocities  values  among  the  solutions.  However, 
the  strong  currents  are  confined  on  a  thin  surface  layer.  Moreover,  SCRUMJ  develops 
sub-scale,  noisy  features  at  the  middle  and  bottom  depth  of  the  water  columns.  Tab.  2. 
summarizes  and  compares  the  model  solutions 

The  across-channel  velocities  also  have  distinct  features.  At  the  vertex  of  the  parabolic¬ 
shaped  shelf,  the  cross-channel  velocities  break  into  in-  and  off-shore  translations  that 
are  associated  with  areas  of  strong  up/downwelling  that  affect  the  whole  water  column 
over  the  continental  slope.  Comparison  of  the  model  solutions  indicates  that  there  is  no 

apparent  correlation  among  the  sub-scale  features  generated  by  each  model. 

•  ^ 

Fig.  2a  compares  the  distribution  of  the  temperature  and  salinity  fields  at  the  middle 
cross  section.  Values  are  approximately  about  7.5  m  below  the  surface.  POM  and  SZM 
temperature  surface  values  exceed  the  maximum  initial  value  over  the  shelf  and  shelf 
break  due  to  the  treatment  of  the  lateral  mixing  operator  (see  Sect.  3.4.).  SCRUMJ 
does  not  present  evidence  of  unrealistic  overheating,  and  the  surface  temperature  is 
the  coolest  among  the  models. 
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4.2.2.  Upwelling 


POM’S  coastal  jet  is  wider  and  located  at  the  mid-continental  slope;  the  mesoscale  ac¬ 
tivity  is  mainly  confined  in  the  upper  water  oolumn  and  generates  at  the  offshore  edge  of 
the  current  system.  At  the  inner  edge,  near-bottom  currents  develop  close  to  the  shelf- 
break,  carrying  bottom  water  over  the  shallow  shelf.  SZM’s  coastal  jet  is  narrower  and 
located  over  the  deepest  part  of  the  continental  slope;  mesoscale  eddies  and  mean¬ 
dering  are  visible  over  the  whole  wind-forced  deep  region.  Similarly  to  the  downweliing 
case,  the  cross-channel  velocity  breaks  into  in/offshore  translations  at  the  shelf  middle 
point  and  no  significant  correlations  are  found  between  the  model  solutions.  The  models 
produce  marked  differences  also  in  the  constituent  fields  and  mixed  layer  configuration. 
Fig.  2b.  compares  the  near-surface  temperature  and  salinity  distributions  at  the  shelf’s 
middle  point.  For  ail  the  models,  near-bottom  currents  at  the  shelf  break  carry  cold  and 
salty  water  into  the  inner  coastal  region.  With  respect  to  the  initial  distribution,  all  models 
present  surface  temperature  minimum  and  salinity  maximum  in  the  proximity  of  the  shelf 
break,  with  a  more  marked  horizontal  gradient  in  the  SZM  and  SCRUM.I  solutions  that 
is  associated  with  a  stronger  local  upwelling  regime.  Over  the  continental  slope,  surface 
temperature  increases  steadly  with  distance  from  the  coast.  ROM’s  initial  surface  val¬ 
ues  and  stratification  are  recovered  outside  the  wind-forced  region.  SCRUM_1  s  surface 
temperature  values  increase  only  sligthly  over  the  deep  part  of  the  channel,  indicating 
that  the  model  is  subjected  to  a  stronger  upwelling  everywhere,  probably  trigged  by  the 
inconsistency  of  the  advective  schemes. 

The  evolution  of  the  salinity  field  is  consistent  with  the  associated  temperature  distri¬ 
bution.  The  less  marked  differences  between  POM  and  SZM  solutions  are  due  to  the 
initial  constituent  distributions  for  which  salinity  is  homogeneous  in  the  upper  30  m,  while 
temperature  is  linearly  stratified. 


4.3.  Internal  Waves  Under  a  Downweliing  regime 

These  experiments  analyze  the  genesis  of  internal  waves  as  a  function  of  the  grid  spac¬ 
ing.  The  problem  is  1-D,  i.e.,  the  bathymetry  is  uniform  in  the  along-shore  direction  and 
ranges  from  20  to  200  m  with  a  mean  slope  gradient  of  0.001 . 

Figs  3a  and  3b  illustrate  the  solution  at  Day  45  for  the  grid  spacing  of  4  and  1  km,  respec¬ 
tively.  Long  runs  are  conducted  to  verify  the  model  stability.  Tests  are  performed  with 


8 


the  old  and  new  version  of  the  SCRUM  model.  It  is  found  that  the  new  version  increases 
dramatically  the  model  solution  and  removes  the  sub-scale  noise  that  characterized  the 
early  program.  Figs  4a-c  illustrate  the  EOF  mode  for  the  solution  at  days  12-15.  All  the 
models  generate  internal  waves  represented  by  the  first  2  modes,  that  account  for  more 
than  50%  of  the  total  variability.  The  remaining  energy  is  distributed  among  the  higher 
mode  differently.  This  is  generally  a  consequence  of  the  different  algorithms  used  by 
the  models  for  the  treatment  of  the  vertical  mixing. 


4.4.  Genesis  and  Evolution  of  Front  Instabilities. 

The  simulations  reproduce  the  test-case  illustrated  in  Wang  (1993,  J.  Mar.  Sci.).  The 
problem  is  configured  on  a  periodic  channei.  The  initial  condition  represents  a  2-D 
density  front  and  the  veiocity  derived  by  each  model  with  diagnostic  calculations.  A 
small  amplitude  pertirbation  is  superimposed  to  the  basic  state.  The  disturbance  starts 
growing  and  the  cyclonic  vorticity  becomes  so  large  that  the  wave  crest  is  bent  backward 
enclosing  an  isolated  eddy  (Fig.  5). 

SCRUM_2  gives  the  more  marked  eddy.  This  is  probably  a  consequence  of  the  constant 
horizontal  eddy  coeficients  (5  m^s— 1)  that  are  equivalent  to  the  mean,  but  also  much 
smaller  than  the  maximum  values  of  the  Smagorinsky  coefficients.  Since  the  Smagorin- 
sky  algorithm  is  a  function  of  the  velocity  shear,  SCRUM_2  is  much  less  diffusive  at  the 
edge  of  the  disturbance.  On  the  other  hand,  SCRUM_2  is  more  diffusive  than  the  other 
models  in  the  lower  region  of  the  channel,  and  the  front  is  smoothed  at  the  warm  side. 

Tests  are  also  performed  to  verify  the  sensitivity  of  the  models  to  the  specification  of  the 
vertical  mixing  algorithms,  and  to  the  vertical  grid  resolution. 


5.  CONCLUSIONS 


This  study  compares  several  ocean  circulation  models  in  coastal  applications.  The  com¬ 
parison  is  restricted  to  3-D,  primitive  equation,  tree-surface  models  on  sigma-like  vertical 
coordinate  systems.  The  models  differ  in  their  treatment  of  the  barotropic  flow,  horizon¬ 
tal  and  vertical  mixing,  temporal  and  spatial  differencing,  and  other  details.  The  physical 
accuracy  and  the  numerical  character  of  the  models  are  evaluated  and  discussed  on  a 
variety  of  test  problems.  The  coastal  study  processes  include  wave  propagation,  up- 
welling  and  downwelling  regimes,  and  genesis  and  evolution  of  fronts. 

One  issue  addressed  is  to  determine  the  robustness  and  portability  of  the  models.  Sim¬ 
ulations  investigate  how  and  if  mixing  algorithms,  filtering  and  smoothing  procedures 
contribute  to  the  models’  numerical  stability  at  the  expenses  of  the  physical  accuracy. 
One  of  the  common  potential  problem  areas  is  the  application  of  the  horizontal  diffu¬ 
sion  operator  acting  on  surfaces  of  constant  sigma.  The  formulation  might  generate 
computational  heating/cooling  sources  in  regions  of  steep  topography. 

The  study  indicates  that  the  solutions  may  sensibly  differ  from  model  to  model.  This 
is  primarily  due  to:  1)  different  response  to  the  calibration  of  the  physical  and  numeri¬ 
cal  parameters,  2)  perturbations  in  the  dynamical  balances.  In  general,  differences  in 
the  primary  features  (such  as  front  location,  strength  of  the  coastal  jet)  can  be  adjusted 
by  adequate  calibrations.  Representation  of  the  secondary  features  (such  as  internal 
wave  propagation,  temperature  and  salinity  anomalies)  may  depend  upon  different  pa- 
rameterizations  and  numerical  treatments  of  the  same  physical  terms  that  may  alter  the 
relative  importance  of  the  dominant  processes. 
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APPENDIX.  The  Asselin  Filter 


Consider  the  advection  equation: 

Ut  =  CUx 

and  the  associated  numerical  approximation: 

U"+l  -  U"-^  ^  (Uj+i  -Uj-i)^  ^ 

2At  2Ax  "" 


Deriving  on  time,  eqs  (A.1)  are  equivalent  to: 


Utt  =c2u. 


Un+l  _  2U“  +  U“-l  _^2  (Uj+i  -  2Uj  +  Uj-ir  ^  ^2  r,n 
At2  Ax2 


The  Asselin  filter  is  applied  over  two  time  steps: 

1)  compute  u“+^  from  eq  (A.1b) 

2)  correct  u®  :  u“  =  u“  +  a(u“’'‘^  —  2u“  +  u““^) 

At  the  following  time  step,  the  variable  u“+2  is  computed  from: 

.nH-2 


U’ 


^  =  C'D”+^  +  0{a) 

2AX  2  /it2 


From  eq  (A.2b)  it  follows: 


Thus,  similarly  to  other  smoothing  procedure,  the  Asselin  filter  is  diffusive  with 
putational  eddy  viscosity  coefficient: 


(A.1a) 

(A.1b) 


(A.2a) 

(A.2b) 


(A.3) 

(A.4) 
a  com- 

{A.5) 
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Tab.  2a.  The  Wind  Driven  Experiment:  Downwelling 
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Tab.  2b.  The  Wind  Driven  Experiment:  Upwelling 
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LATERAL  MIXING 
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Fig.  1.  Cross  channel  sea  surface  temperature  from  the  a-level  and  geopotential  horizontal 
mixing  algorithms  (a).  Sea  surface  temperature  differences  between  the  two  solu¬ 
tion  (b).  Along-shore  coastal  jet  structure  for  the  two  lateral  mixing  formulations  {c 
and  d). 
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